
% load BHMR_RESTAT_STRUCTURAL.mat
% 
 Popsize = 1729
 mFixed=15
% 
% % treatment = 1  --> Discount only
% % treatment = 2  --> Fixed commitment, no discount
% % treatment = 3  --> Fixed commitment with discount
% % treatment = 4  --> Flexible commitment, no discount
% % treatment = 5  --> Flexible commitment with discount
% % treatment = 6  --> Control
% 
% TreatmentPop = NaN(Popsize,1);
% 
% for i = 1:Popsize
% 	if treatment(i,1) == 6    
% 		TreatmentPop(i,1)= 0;
%     else
%         TreatmentPop(i,1)=treatment(i,1);
%     end
% end
% 
% % TreatmentPop = 0  --> Control
% % TreatmentPop = 1  --> Discount only
% % TreatmentPop = 2  --> Fixed commitment, no discount
% % TreatmentPop = 3  --> Fixed commitment with discount
% % TreatmentPop = 4  --> Flexible commitment, no discount
% % TreatmentPop = 5  --> Flexible commitment with discount
% 
% 
% attendsCamp = anyvisit;
% 
% %Should this be 0 or NaN for control and discount only?
% acceptsContract = zeros(Popsize,1);
% 
% for i = 1:Popsize
% 	if TreatmentPop(i,1)~=0 & TreatmentPop(i,1)~=1 
%         acceptsContract(i,1)=takeup_v2(i,1);
%     end	
%     if acceptsContract(i,1)==Inf
%        acceptsContract(i,1)=NaN;
%     end
% end
% 
% commitmentAmountChosen = NaN(Popsize,1);
% 
% 
% for i = 1:Popsize
% 
% 	if (TreatmentPop(i,1)==4 | TreatmentPop(i,1)==5) & acceptsContract(i,1)==1
%         commitmentAmountChosen(i,1)=ccamount(i,1);
%     end
%     if (TreatmentPop(i,1)==2 | TreatmentPop(i,1)==3) & acceptsContract(i,1)==1
%         commitmentAmountChosen(i,1)=mFixed;
%     end	
%      if commitmentAmountChosen(i,1)==Inf
%         commitmentAmountChosen(i,1)=NaN;
%     end
% end
% 
% 
% Distance = NaN(Popsize,1);
% 
% for i = 1:Popsize
% 
% 	if distance_ehp(i,1)==1
%         Distance(i,1)=.5;
%     end
%     if distance_ehp(i,1)==2
%         Distance(i,1)=3;
%     end
%     if distance_ehp(i,1)==3
%         Distance(i,1)=6;
%     end
% 		
% end
% 
% %recode vars needed for heterogeneity inf to nan
% for i = 1:Popsize
% 
% 	if chronic_dis(i,1)==Inf
%         chronic_dis(i,1)=NaN;
%     end
%     
%     if times_sick_q4(i,1)==Inf
%         times_sick_q4(i,1)=NaN;
%     end
%     if hibp(i,1)==Inf
%         hibp(i,1)=NaN;
%     end
%     if s2_hyper(i,1)==Inf
%         s2_hyper(i,1)=NaN;
%     end
%     if s1_2_hyper(i,1)==Inf
%         s1_2_hyper(i,1)=NaN;
%     end
% 	if hh_income(i,1)==Inf
%         hh_income(i,1)=NaN;
%     end	
%     
%     if taking_meds(i,1)==Inf
%         taking_meds(i,1)=NaN;
%     end	
%     if gender(i,1)==Inf
%         gender(i,1)=NaN;
%     end
%     if employee(i,1)==Inf
%         employee(i,1)=NaN;
%     end
%     
%    if pbias(i,1)==Inf
%         pbias(i,1)=NaN;
%    end
%     if sickv(i,1)==Inf
%         sickv(i,1)=NaN;
%     end
%     if literate(i,1)==Inf
%         literate(i,1)=NaN;
%     end
%     if trusts_ehp(i,1)==Inf
%         trusts_ehp(i,1)=NaN;
%     end
%      if hbp_freq_check(i,1)==Inf
%         hbp_freq_check(i,1)=NaN;
%      end
%      if procrastinate(i,1)==Inf
%         procrastinate(i,1)=NaN;
%      end
%      if impatient(i,1)==Inf
%         impatient(i,1)=NaN;
%     end
% end
% %%% Make sure hbp indicators match up
% for i = 1:Popsize
% 
%     if s2_hyper(i,1)==1
%         hibp(i,1)=1;
%     end
%     
%     if s1_2_hyper(i,1)==1
%         hibp(i,1)=1;
%     end
%     		
% end
% %%% Generate indicators
% sickness_ind=chronic_dis+times_sick_q4;
% hbp_ind=hibp+s1_2_hyper+s2_hyper;
% trust_ind=trusts_ehp+literate+sickv;
% aware_ind=procrastinate+impatient;
% 
% 
% save realdata_new.mat TreatmentPop acceptsContract commitmentAmountChosen attendsCamp Popsize Distance sickness_ind hbp_ind hh_income taking_meds employee pbias trust_ind hbp_freq_check aware_ind gender literate

% %% Estimation
% load realdata_new.mat
% 
Popsize = 1729;
% 
% %%% Now set fixed parameters for model / estimation
% 
 f	=	30;          % Per visit monetary fee to be paid 
 d	=	15;          % Per visit discount offered in some treatments - 
 mFixed   =  15;     % Per visit commitment amount -- Amount person gets back at time of visit (since fee paid up front under all commitment contracts) 
 mFlexMin = 0;        % The minimum commitment we required in the flexible commitment amount (money you get back at time of service for showing up)
 Sim = 50;
 cSD=0;
% 
% %%% Now, set standard draws from two normal distribution to hold fixed throughout parameter search for random coefficients
% 
 mu = [0,0,0,0];
 sigma = eye(4);
 RC = mvnrnd(mu, sigma, Sim);
 J = unifrnd(0,1,1,Sim);
% 
% %%% Have everything we need now to do constrained non-linear optimization with simulated maximum likelihood
% 
% %Heterogeneity vars and indicators
 DistanceData = Distance*ones(1,Sim);
 sickness_indData = sickness_ind*ones(1,Sim);
 hbp_indData = hbp_ind*ones(1,Sim);
 hh_incomeData = hh_income*ones(1,Sim);
 taking_medsData = taking_meds*ones(1,Sim);
 employeeData = employee*ones(1,Sim);
 literateData = literate*ones(1,Sim);
 pbiasData = pbias*ones(1,Sim);
% 
% 
 trust_indData = trust_ind*ones(1,Sim);
 hbp_freq_checkData = hbp_freq_check*ones(1,Sim);
 aware_indData = aware_ind*ones(1,Sim);
 genderData = gender*ones(1,Sim);
% Contract and attendance vars

% attendsCampData = attendsCamp*ones(1,Sim); 
% commitmentAmountData = commitmentAmountChosen*ones(1,Sim);
% acceptsContractData = acceptsContract*ones(1,Sim);

%% PRIMARY ESTIMATION ROUTINE

%  parameters0 = [0.7 0.5 5 5 0.5 10 5 5 5 0 0 5 5 5];   % Starting points for global search of the params being optimized over
%  lb = [0.001 0.001 0.001 0.001 0.001 0.001 -50 -50 -50 -0.5 -0.5 -50 -50 -50];  % lower bound of the params being optimized over
%  ub = [1 2 100 100 1 200 50 50 50 0.5 0.5 50 50 50];  % upper bound
%  gs = GlobalSearch;  % global search problems search cleverly over the parameter space with multiple starting points
%  options = optimset('MaxFunEvals',40000,'MaxIter',40000,'TolFun',10^-7);
%  likelihood =  @(parameters) Likelihood_woDistance([parameters(1) parameters(2) parameters(3)  parameters(4) parameters(5) 30 0 parameters(6) 30 parameters(7) parameters(8) parameters(9) parameters(10) parameters(11) parameters(12) parameters(13) parameters(14)],TreatmentPop, RC, J, Popsize,f,d,mFixed,mFlexMin,Sim,acceptsContractData,commitmentAmountData,attendsCampData,hbp_indData,sickness_indData, employeeData, taking_medsData, genderData, literateData);   % note that some params are "fixed" to exogenous values
%  problem = createOptimProblem('fmincon','x0',parameters0,'objective',likelihood,'lb',lb,'ub',ub,'options',options);
%  [xmin,fmin,flag,outpt,allmins] = run(gs,problem);

 

